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Abstract 

We present a detailed study of a 4 parameter family of elliptic weights on tilings of a hexagon 
introduced by Borodin, Gorin and Rains, and generalize some of their results. In the process, we 
I connect the combinatorics of the model with the theory of elliptic special functions. 

We first analyze some properties of the measure and introduce canonical coordinates that are 
I I useful for combinatorially interpreting results. We then show how the computed n-point function 

i-S^ (called the elliptic Selberg density) and transitional probabilities connect to the theory of BC„- 

symmetric multivariate elliptic special functions and difference operators discovered by Rains. In 
fH particular, the difference operators intrinsically capture the combinatorial model under study, while 

4-^ the elliptic Selberg density is a generalization (deformation) of probability distributions pervasive in 

the theory of random matrices and interacting particle systems. 
^ Based on quasi-commutation relations between elliptic difference operators, we construct certain 

' ' natural measure-preserving Markov chains on such tilings. We then immediately obtain and describe 

I an exact sampling algorithm from such distributions. We present sample random tilings from these 

^ measures showing an arctic boundary phenomenon. Interesting examples include a 1 parameter 

family of tilings where the arctic curve acquires 3 nodes. 

Finally, we show that the particle process associated to such tilings is determinantal with cor- 
^-H relation kernel given in terms of the univariate elliptic biorthogonal functions of Spiridonov and 

Zhedanov. 
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1 Introduction 

This paper examines work began by Borodin, Gorin and Rains in |BGR10| . In op. cit., the authors 
examined q-distributed boxed plane partitions from several perspectives, but the g-distributions were 
obtained as limits of the elliptic distribution briefly appearing in the Appendix. The present paper takes 
the Appendix of |BGR10| and expands upon it, following the steps in B GRIO] and |BG09] . However, 
since we are working at the elliptic level (rather than a degeneration as in iBGR10| ). new tools are needed 
to generalize the results of [BGRlOj . These tools belong to the area of elliptic special functions, an active 
area of research in algebra and analysis generalizing, among other things, the Askey and g-Askey schemes 
of orthogonal polynomials (as described in [KS, for example). Thus, in some complementary sense, while 
being a generalization of [ BGRlOj , the paper is an application of multivariate tools introduced by Rains 
in |RailO| and jRai06j (the first is more analytic, the second being more algebraic) which build upon 
univariate elliptic biorthogonal functions found by Spiridonov and Zhedanov a few years earlier in |SZOO| . 
Work in the area of elliptic special functions started with Frenkel and Turaev's discovery of elliptic (theta) 
hypergeometric series ( |FT97) ) - the authors of op. cit. cite Baxter's work (see his book |Bax82] ) as the 
genesis of the theory. 

The history of the problem starts with random uniformly distributed boxed plane partitions. Much 
is known about these: asymptotics and frozen boundary behavior f |CKP01) . |CLP98| . |KU07j ): cor- 
relation kernel via Hahn orthogonal polynomials (see |Joh05) , [BGOQj . |Gor08p : exact sampling algo- 
rithms ( |BG09| V Somewhat central to the subject is the topic of discrete Hahn orthogonal polynomials 
(which themselves are terminating generalized hypergeometric series). One level up and we arrive at q- 
distributions on boxed plane partitions ( [BGRlOj . and |KO07| for the variational problem used to derive 
the limit shape for the qri^o''""^ distributions). Almost as much as above is known about these, and 
central to the subject are certain discrete q-orthogonal polynomials (g-Racah, g-Hahn) from the g- Askey 
scheme, which themselves are terminating g-hypergeometric series (see jGR04j for a full description, [KSj 
for a distillation of the formulas). 

The present work analyzes the elliptic level (the distribution was introduced in the Appendix of 
[BGRIO] . but also independently from a slightly different perspective in [Sch07] ). We look at two 
aspects: exact sampling algorithm and correlation kernel. The third aspect in |BG09j and [BGRlOj is 
obtaining asymptotics of the correlation kernel and through this obtaining the frozen boundary behavior 
in the large scale limit. While we indeed see a frozen boundary behavior in our case and can characterize 
it via variational techniques (and we present computer simulations of the results), we cannot yet analyze 
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the asymptotics of elliptic biorthogonal functions (techniques used in previous works - e.g., in [BGRIO] - 
fail if we replace orthogonal polynomials by elliptic biorthogonal functions) . More direct techniques like 
solving the variational problem described in |KQ07| for the g-Hahn case and in Section 2.4 of [BGRIO] 
seem computationally intractable so far. The reason is the associated complex Burgers equation one has 
to solve becomes considerably more complicated. Nevertheless, it is a (new) feature of the elliptic model 
that the frozen boundary can have 3 nodal points (as seen in the computer simulations). 

From a different perspective, we try to create a bridge between elliptic special functions discussed in 
the references given above and combinatorics of tilings of hexagons (equivalently, dimer coverings of the 
appropriate graph etc. See Section 2 for interpretations). To wit, we give a combinatorial interpretation 
to several objects appearing in the theory of elliptic special functions: the (t = q case) multivariate 
elliptic difference operators discovered by Rains ( IRailOj ). the A-symbols of |Rai06| and the (univariate) 
elliptic biorthogonal functions of Spiridonov and Zhedanov f |SZOO| ). 

This paper tries to emulate the organization of |BG09| and [BGRlOj . but with notation heavily 
influenced by |Rai06] . It is organized as follows: in the remainder of the introduction, we set up most 
the important notation and terminology. 

We set up the combinatorial and probabilistic aspects in Section 2 (in the Appendix we consider a 
different way of assigning weights to rhombi that is manifestly more symmetric. We prefer to use the 
formulas from Section 2 though, at the cost of symmetry breaking since they lead to shorter computations 
and arguments). Also in this section we study positivity of our a priori complex measure and introduce 
various coordinate systems used throughout the paper, including the important canonical coordinates 
(which embed our model in a certain square of an elliptic curve). 

In Section 3 we compute relevant distributions and transition probabilities. Sections 2 and 3 are an 
expansion and in depth analysis of the Appendix in [BGRIO] . 

Section 4 recalls some definitions and properties of elliptic tools introduced by Rains ( [RailO[ . [Rai06] 
- we refer the reader to these works for the proofs we omit) and then connects these with the probability 
and combinatorics being studied. We show that the constraints of the model are intrinsically captured 
by the elliptic difference operators under discussion. 

Section 5 describes a perfect sampling algorithm for such elliptic distributed boxed plane partitions. It 
is based on the idea of forming a new measure-preserving Markov chain out of two old quasi-commuting 
ones (as in [BFlOj : see also [DF90j ). The algorithm starts from a deterministic parallelogram shape 
and samples relatively easy distributions to successively transform the parallelogram into a hexagon 
accordingly distributed (by increasing one side by 1, and decreasing another side by 1; a parallelogram 
can be seen as a hexagon with two sides of length 0). We use the quasi-commutation relations for the 
elliptic difference operators of Section 4 to construct this algorithm. 

Section 6 deals with the correlation kernel. We start by recalling facts about univariate elliptic 
biorthogonal functions and show that the time increasing (decreasing) Markov process is indeed deter- 
minantal, with correlation kernel given as a determinant of (a matrix composed of) elliptic biorthogonal 
functions. These replace the orthogonal polynomials discussed above. 

In Section 7 we present some computer simulations obtained from the algorithm described in Section 
5. Choice of parameters for obtaining the trinodal cases (surfaces where the arctic circle has 3 nodes at 
3 vertices of the hexagon) are also explained. 

We end with the Appendix, which provides a highly symmetric view of the entire picture 

For the remainder of the section, we will set the notation that will appear in the rest of the paper. 
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We define the theta and elUptic Gamma fmictions f lRuiQ?] ) as follows: 

A:>0 



. n ^ ' ^ 



fc,/>o ^ ^ 

Note the elliptic gamma function is symmetric in p and q. The theta-Pochhammer symbol (a gener- 
alization of the g-Pochhammer symbol) is defined, for m > 0, as 

0<i<m 

As is usual in this area, presence of multiple arguments before the semicolon (inside theta or elliptic 
Gamma functions) will mean multiplication. To wit: 

6'p(uz='=^; g),„ = 0p{uz; q)mOp(u/z; q)^; T^p,q{a, b) = Vp^q{a)Tp^q{b). 

We have the following important identities {n > an integer): 

dp{x) = 0p{p/x) 

Op{px) = ep{l/x) = -{l/x)9p{x) (1) 
Tp^q{q'^x) = 9p{x;q)„Tp^q{x) 

The last identity in ([T]) can be extended for rt < or even for non integer n to provide a generalization 
of the theta-Pochhammer symbol for negative or even non-integer lengths. These identities also extend 
to the following among theta-Pochhammer symbols: 

9p{a; q)n+k = Op{a; q)nOp{aq"';q)k 
ep{a;qU^0p{q^-"/a;qU^arqi"-) 

dp{a,q)n-k - ^ I ^( ) (T^^' 

0p(q /a;q)k a 

n „^ ^p{a\(l)kOp{qla\q)n ^-nk (2) 

'^^"'^ '^^'== ep{q^->^la-q)^ 

Op[aq ";(?)„ Op{q/a;q)n a 

a f n 9p{a;q)k0p{aq'';q)n dp{a;q)n+k 
9p(aq -qjk = r = r — 

We will use the above identities throughout for simplifying computations without explicitly referring to 
them. 
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If /(a;i, Xn) is a function of n variables defined on (C*)", we call it BCn- symmetric if it is symmetric 
(does not change under permutation of the variables) and invariant under Xk — >■ l/a;^ for all k. We will 
call it a BCn-symmetric theta function of degree m if in addition, it satisfies the following: 

fipxi, ...,Xn) = (^)™/(a;i, ■■•,a;„). 
pxf 

The prototypical example of a _BC„-symmetric theta function of degree 1 is: 

n 

l<k<n 

We now define the following function (which will play an important subsequent role): 

ip{z,w) = z^^9p{zw, z/w) (3) 

Note (fi is _BC2-antisymmetric {ip(z,w) — —ip{w,x)) of degree 1. It is a consequence of the addition 
formula for Riemann theta functions that 



'P{z,x) _ f{z,y) \ if{w,x)ip{w,y) 
(p{w,x) (p{w,y)J ip{z,w) 



for arbitrary z, w. We observe that the expression in parentheses appearing above is a Vandermonde-like 

Y - f'^''''^^ so uy(7- 7-') 



factor in transcendental coordinates X = ^4^^,y = ^t^^, so if{zi,Zj) is an "elliptic analogue" of the 



( Vandermonde) difference Zi — zj . This is indeed the case if one takes the right limit and then a product 
over i < j. To wit: 

limp^o ¥'(g''' , g^O 

Imi — Xi — Xj. 

Notationally, for a function / of n variables, we will use the abbreviation f(...Xk---) to stand for 

f{xi , . .., Xji) ■ 

We will make reference to the delta symbols defined in |RailO| (see also |Rai06| - we are in the case 
t = q in the notation from both references), which we define here. We fix A € to" a partition (that is, a 
partition with at most n parts all bounded by to). Define the partition 2A^ by (2A^)i — 2(A|-i/2]). 



elix;q)^Y[epiq'-'x;q)x, 

i<i 

l<i 

p+( N TT dp{q^-'-^x;q)x,+x, yr Op{q^-'-^x) -pr 6lp(g^-^'x; g)2A, 
' i4 y ^^,(g2— .+A.+A,^) 11 ep{q^- -^-^x- q),^ 

eA-(-;'^)= n yf:::^'^^;-^- -U o^'^^^l , n^.(^"-'-;^)^. 



AA(a|...6i...;g) 



e°{...h..-q) e°2X-{pqa;q) 



eO(...?^...;q) e^(pq,<j;g)e+(pa,a;g) 
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Of interest will be the A-symbol with six parameters to, ^i, ^2, ^3, uo, ui satisfying the balancing 
condition (?^"~^toii^2t3'«oUi — q- Because the usual balancing condition has pq on the right hand side 
(the reader should consult the Appendix of |RailOj for more on why this is necessary), we multiply ui 
by p (this choice is arbitrary, so a priori some symmetry is broken, but this will not affect our results) . 
We define the discrete elliptic Selberg density as: 



A^{q^^-hl\q^,q^-hoti,q^-hot2,q''-Hots,q''-%Uo, g"-4o(pui); q) = const ■ l[{ip{z,, z,))^ (4) 



nh(2n~l)n ( 2\ ^p(^0: ^O^li ^0^2, ^0^3, tpUQ^tQUi] q)l^ 

where li ^ n — i + Xi, Zi — q^'^to and the constant is independent of A and present in the formula to 
make the A-symbol elliptic in all of its arguments (the value of the constant can be computed, but such 
constants will be immaterial for the rest of the paper) . This discrete elliptic Selberg density is the weight 
function for the discrete elliptic multivariate biorthogonal functions defined in |Rai06| . Notice it can be 
written more symmetrically in terms of the Zi 's and the elliptic Gamma functions as 

2n-l/) / 2\ ^p,q{toZi7tlZi7t2Zi,t3Zi,UoZi,UiZi) 



const ■ Y[i^iz,, zj)f -Xlzf'-^epiz. 



r {-^Z- ^Z- -3-z- -2-z-V 



i<j i "P^^Vto"*' t2-»' t3"«' «0 "1 

We will denote by E the elliptic (Tate) curve C*/{p) for some \p\ < 1. An elliptic function / (of 1 
variable) will just be a function defined on E (that is, f{px) = f{x)). 

Throughout the remainder, constants (by which we mean factors independent of the variables usually 
denoted by Xk,yk, Zk) will largely be ignored (and we will write const wherever this appears), but they 
are there to make measures into probability measures (i.e., normalizing factors) or to make certain 
functions elliptic (i.e., invariant under p-shifts). Their values can often be recovered, and we comment 
on how to recover them whenever possible. 

Finally, throughout this paper we will freely use two different systems of coordinates for our model 
(related by a simple affine transformation - see the next section). While this may seem redundant, coor- 
dinatizing in two different ways will more aptly reveal different features of the elliptic special functions 
and difference operators under study. 

Acknowledgements 
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2 The model 
2.1 Interpretations 

We consider random tilings of an a x x c regular hexagon embedded in the triangular lattice (with 
Cartesian coordinates by tiles of three types, as can be seen in the Figure [T] The probabilistic 
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details are set out in Section 2.2 We will find it more convenient to encode the hexagon via the following 
three numbers: 

N ^ a,T ^b + c,S ^ c. 




Figure 1: A tiling of a 3 x 2 x 3 hexagon and the associated stepped surface 

Equivalently, these tilings can be thought of as dimer matchings on the dual honeycomb lattice (every 
rhombus in a tiling is a line matching two vertices in the dual lattice), stepped surfaces, boxed plane 
partitions (6 x c rectangles with positive integers < a filled in that decrease weakly along rows and 
columns starting from the top left corner box) or 3D Young diagrams (any section parallel to one of the 
three bounding walls is a Young diagram). 

A yet different way of viewing such tilings, important hereinafter, is as collections of non-intersecting 
paths in the square lattice. The paths start at N consecutive points on the vertical axis (counting from 
the origin upwards) and end at N consecutive points on the vertical line with coordinate T. Each path is 
composed of horizontal segments or diagonal (Southwest to Northeast, slope 1) segments, and the paths 
are required not to intersect. Figure [2] explains this, and also introduces the coordinate frame (t, x) that 
will be used for computational convenience in various sections to follow: 

= {t,x~t/2). 

Following the notation in ;BG09], let n{N,S,T) denote the set of N non-intersecting paths in the 
lattice starting from positions (0, 0), (0, - 1) and ending at positions (T, S), {T,S + N - 1). 
Each path has segments of slope or 1 (paths go either horizontally or diagonally upwards from left to 
right). Set 

X^_*j, ^{xeZ: max(0, t + S - T) < x < mm{t + N-1,S + N-1)} 
X^'lp = {X = {xi,...,xn) e (X;^*j,)^ : xi < X2 < ... < xn}. 

S t • 

rp is the set of all possible particle positions in a section vertical section of our hexagon with horizontal 

coordinate t (in (i, x) coordinates). X^'^ is the set of all possible A^-tuples of particles in the same vertical 
section. 

For X G n{N, S, T), we have X = {X{t))o<t<T and each X{t) e X^'^. X is a discrete time Markov 
chain as it will be shown. 
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Figure 2: Duality between tilings and non-intersecting paths 
2.2 Probabilistic model 

We will now define the probability measure on il.{N, S, T) that will be the object of study. For a tiling 
T corresponding to an X e ^{N, S, T) we define its weight to be: 

w{7) = J] w{l) 

^tz{ horizontal lozenges} 

where by a horizontal lozenge we mean a lozenge whose diagonals are parallel to the i and j axes. The 
probability of such a tiling would then simply be: 



Prom - ^ "^^^ 
The weight function w on horizontal lozenges I is defined by 



w{l) 



where is the coordinate of the top vertex of the horizontal lozenge I, ui,U2,q,p are complex 

parameters {\p\ < 1) and ui — q~^vi, U2 = ^2 (the reason for this break in symmetry is that it will make 
other formulas throughout the paper more symmetric). 

Remark 2.1. Only considering weights of horizontal lozenges for a tiling of a hexagon is equivalent to 
considering all types of lozenges but assigning the other two types weight 1 (i.e., each lozenge that is 
not horizontal has weight 1). This is a break in symmetry that can easily be fixed. However, for the 
remainder of the paper we prefer this non-symmetric weight assignment system as it makes computations 
easier. Nevertheless, we show in Appendix [8] that we can assign weights to the 3 types of lozenges in 
a S'3-invariant way (i.e., invariant under permuting the 3 types of lozenges or equivalently the 3 spatial 
directions). 
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Figure 3: Going from 3 dimensions to 2 dimensions 




Figure 4: A full 1x1x1 box (left) and an empty one (right) 



This weight on dimer coverings of a hexagon was derived in |BGR10| (see also |Sch07| for elliptic 
enumeration of lattice paths) . 

The connection with elliptic functions will now be explained. Fix a horizontal coordinate i, denote 
by w[i,j) the weight of the horizontal lozenge with top vertex coordinates (j,j), and observe that for 
two consecutive vertical positions we have (1/1^2^3 ~ 1): 

W{i,j-1) e'p(9J-3i/2+ly^^qj+3»/2+ly2,5-2i+ly3) 
_ g3gp(gJ-3V2-S-l^^^^j-+3.:/2~l^^^g-2j-+5-l/^^^^) 

~ 6lp(gJ-3V2-S+l„^^qJ+3^/2+l„2,g-2J + S+l/^;^„2) 

In 3-dimensional coordinates (x, y, z) pictured in Figure [3] (note we only consider surfaces in 3 di- 
mensions that are stepped, meaning there is a 1-1 correspondence between the 2D tiling picture and the 
3D surface picture) with i = x — y,j = z — {x + y)/2, the weight ratio looks like 



^ ^ ^ w(full box) ^ q^ep{ui/q,U2/q,U3/q) 
w(empty box) 6p{uiq,U2q,U3q) 



where 



and {x, y, z) is the 3-dimensional centroid of the 1x1x1 full cube (on the left in Figure |4]) with top lid 
the horizontal lozenge with top vertex coordinate («, j). 
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The word elliptic now becomes clear as r in ([8]) is an elliptic function of q (that is, defined on E - 
see the Introduction for details). Moreover, r is the unique elliptic function of q with zeros at ui,U2,U3 
and poles at 1/ui, 1/m2, l/wa normalized such that r(l) = 1. Of interest is also that r is elliptic in Uk 
for fc = 1, 2, 3 subject to the condition that nfe=i — 1- 

Remark 2.2. r is invariant under the natural action of 6*3 permuting the u^'s (and of course the 3 axes: 
x,y,z). 

Wc can view our tilings as stepped surfaces composed of 1 x 1 x 1 cubes bounded by the 6 planes 
X — 0,y = 0, z = 0,x = b,y = c, z = a. Then the 2 dimensional picture in Figure [T] can be viewed as a 
projection of the 3 dimensional stepped surface onto the plane x + y + z — 0. 

For T a tiling, we have 



wt{7)= Yl w{i,j) 

G T 

where (i,j) are the coordinates of the top vertex of a horizontal lozenge. Grouping all 1 x 1 x 1 cubes 
into columns in the z direction with fixed {x^y) coordinates (see Figure |3|, we obtain: 

wt{7) = const ■ n ""^''-^^ , 
w(z,j - 1) 

where the product is taken over all cubes (visible and hidden) of the boxed plane partition and («, j) is 
the top coordinate of the bounding hexagon of a 1 x 1 x 1 cube. Note to get to this equality we have 
merely observed that wt(empty box) is a constant independent of i and j. We can further refine this 
(rearranging the terms in the product and gauging away more constants - see Section 2.3 of [BGRlOj 
for more details) as: 

wt{7) = const ■ ^ ^"^^Y^i^ ' ^ = const • nr(i,j)''^"^ 

where v = {xo,yo,Zo) ranges over all vertices on the border (but not on the bounding hexagon) of the 
stepped surface with xo,yo,ZQ integers (equivalently, v ranges over all vertices of the triangular lattice 
inside the hexagon, but we view w in 3 dimensions). h{v) is the distance from v to the plane x + y + z = 
divided by a/3 : h{v) — {xn + yo + zo)/3. 

2.3 Positivity of the weight 

The content of the previous subsection shows that in order to make the whole model well defined as a 
probabilistic model, it suffices to establish positivity of the elliptic weight ratio r(i,j) — w{i, j)/w{i, j — 1) 
defined in ([?]) (where (z, j) is the location of a given horizontal tiling and ranges over all possible horizontal 
tilings inside the hexagon). Recall that 

/. q^0p{ui/q,U2/q,U3/q) 
9p[qui,qu2,qu3) 
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where iii — (/■'"^'/^mi, M2 = g''+^*/^W2, "3 — q~'^^uz and U1U2U3 = 1. We recall that r is elliptic in Uk for 
k = 1, 2, 3 as well as in q. In order to make r positive, we will first restrict ourselves to the case where r 
is real valued. This means r is defined over a real elliptic curve, and we have — l<P7^0<l(a priori, 
p is complex of modulus less than 1; p G (—1, 1) — {0} is equivalent to E being defined over M - for more 
on real elliptic curves, we refer the reader to Chapter 5 of |Sil94| ). We then ensure positivity of r by an 
explicit computation. We will of course have two cases: p < and p > 0. We deal with the case p > 
throughout (and make remarks when necessary for p < 0). 

Now that we have restricted ourselves to real elliptic curves E, we first note that q G E (i.e., r is 
elliptic as a function of q). For a chosen < p < 1 there are two non-isomorphic elliptic curves defined 
over K (since Gal{C/M.) = Z/2Z), both homeomorphic to a disjoint union of two circles (every real 
elliptic curve is topologically homeomorphic to a circle if p < or with a disjoint union of two circles if 
p > - one can just see this by plotting the Weierstrass equation in and compactifying) : 

E =R R*/p^ and 

E=R {ueC7p^:|7/pe{l,p}} 

We will call the first case real and the second trigonometric (abusing terminology, since both are real 
elliptic curves). We will analyze the trigonometric case, but the real case is similar. In the trigonometric 
case, the curve has two connected components (circles): the identity component (it contains the points 
1 and -1) and another component that contains the other 2-torsion points: ^^/p■ There will be 3 cases 
to be analyzed which we list now and motivate after (if p < there is only one component so the 3 cases 
coalesce to only one - Case 2.): 

• Case 1. q lies on the non identity component {\q\ — ^/p). 

• Case 2. q and all the u^'s (and so the Wfc's) lie on the identity component {\q\ — \ui\ ~ \u2\ = 

= 1) 

• Case 3. q and one of the u^'s lies on the identity component, the other two u^'s lie on the non 
identity component 

To analyze positivity at a fixed site («, j) inside the hexagon, we note that r(g) has zeros at and 
poles at 1/uk {k = 1,2,3). We note r = ±1 at g = ±1 so at least one Uk (along with its recipro- 
cal/complex conjugate 1/uk) needs to be on the identity component (so that r can change signs on the 
identity component). Since r = —1 at q = and U1U2U3 = 1, either exactly one or all three of the 

w's need to be on the identity component. This motivates the three choices above. 

Case 1. will never lead to positivity for all four admissible sites inside a 1 x 2 x 2 hexagon (see 
Figure [s]) , so we can eliminate it (if a 1 x 2 x 2 hexagon is never positive, much larger ones which are of 
interest to us will also never be as they contain the 1x2x2 case). For a proof, we suppose that ui is 
on the identity component, and 1*2, M3 are (along with q) on the non identity component (the case where 
all three u's are on the identity component is handled similarly). The u's differ from the u's by integer 
powers of q given in the last three columns of the following table (for the four admissible (i, j) pairs in 
the 1x2x2 hexagon): 



j 


i 


J - 3^/2 J - 


h3i/2 


-2j 


1/2 


1 


-1 


2 


-1 


1 


2 


-2 


4 


-2 





2 


-3 


3 





1/2 


3 


-4 


5 


-1 
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Figure 5: The admissible sites inside a 1 x 2 x 2 hexagon 



Notice mod 2 (and we only care about mod 2 as is on the identity component), the four vectors 
(from the last 3 columns of the table) above are (1, 0, 1), (0, 0, 0), (1, 1, 0), (0, 1, 1). The corresponding 
Mfc's we get by multiplying each Uk by q to the power coming from the vector (0, 1, 1) - (■ui,{(2, W3) = 
(g~^wi, (7^W2, (/"^wa) - will all be on the identity component, which means the elliptic weight ratio will 
be negative at the site (j,j) = (1/2,3) as q is on the non identity component. This is a contradiction. 
The other cases are handled similarly, leading to contradictions. This proves q must be on the identity 
component, so only cases 2. and 3. above can lead to positive hexagons. 




Figure 6: The identity component of E =r {u € C*/p^ : |up G {l,p}}. For positivity of r throughout 
the hexagon (i.e., for all admissible Wfc's), q must always be closer to 1 than any u^^ as depicted. 

I will next discuss the case where q and all are on the identity component (case 2. above; for case 
3. the reasoning is similar). For a fixed site inside the hexagon, the 3 ii^s and their reciprocals 
(complex conjugates) break down the unit circle into 6 arcs (see Figure |6]) and q must be on one of the 
three arcs where r is positive (as depicted in the figure). If we want to ensure positivity of the ratio for 
all 4 admissible sites (z, j) within a given 1x2x2 hexagon (Figure [5|, we first observe that for \x\ = 1: 

0,{x) = {l-x)\{\l~p^x\\ 

i>l 

So we reduce to positivity of the corresponding four functions Jli=i 2 3 ^I'-r/q ■ Through standard trigono- 
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metric manipulations wc thus want positivity of each of the following functions: 

sin7r(Q;i — a) sin7r(Q;2 — a) sin7r(a3 — a) 

sin7r(Q;i + a) sin7r(Q;2 + a) sin7r(a3 + a) 

sin7r(ai) sin7r(a2 — 3a) sin7r(a3) 

sin7r(ai + 2a) sin7r(Q;2 — a) sin7r(a3 + 2a) 

sin7r(Q;i — 3a) sin 77(02) sin7r(a3 + a) 

sin7r(Q!i — a) sin7r(Q!2 + 2q:) sin7r(Q:3 + 3q;) 

sin7r(Q;i + a) sin7r(Q!2 — 2q:) sin7r(Q:3 — 2a) 

sin 7r(Q;i + 3a) sin7r(a2) sin7r(a3) 

where 27rai = argUi,ai +02+03 G {0, 1, 2}, 27ra = argq and (a, oi, 02) are defined on the 3-dimensional 
unit torus K'^/Z'^. If we restrict to the fundamental domain [0, 1]'^ and look at all the regions (polytopes) 
cut out by the planes (linear functions) in the arguments of the sines above (divided by tt), we find (by 
solving the appropriate linear programs via Mathematica) that there exists only one region of positivity 
for all 4 functions. We can characterize the region best in terms of Figure |6] That is, as (i, j) range over 
all 4 sites inside a 1 x 2 x 2 hexagon, there should not be any Uk (k = 1,2,3) or any ii^^ on the arc 
subtended by 1 and q (and that does not contain -1). 

Remark 2.3. In view of the above, for any choice of a reasonably large hexagon (say one that contains 
a 1 X 2 X 2 hexagon) and parameters Mi, 1*2,^3 (satisfying the balancing condition), the set of q giving 
rise to nonnegative weights is a symmetric closed arc containing 1. 

2.4 Degenerations of the weight 

Certain degenerations of the weight have been studied before (among the relevant sources for our purposes 
are |BG09| . jBGRlOj . |Joh05| . |KO07| . |Gor08p from many angles. For example, when g = 1 the weight 
in ([6| becomes a constant independent of the position of the horizontal lozenges, and so we are looking 
at uniformly distributed tilings of the appropriate hexagon. An exact sampling algorithm to sample such 
a tiling was constructed in |BG09| and the theory behind this (as well as behind other results for such 
tilings) is closely connected to the theory of discrete Hahn orthogonal polynomials (see |Joh05| , |BG09) , 
[Gor08| l. The frozen boundary phenomenon (the shape of a "typical boxed plane partition") was first 
proven in |CLP98| and then via alternate techniques (and generalized) in jCKPOl] and |KO07j . 

A more general limit than the above is the following: in ([6]) we let vi = V2 = Hy/p and then let 
p — > 0. This is the g-Racah limit (named after the discrete orthogonal polynomials that appear in the 
analysis). This limit is the most general limit that can be analyzed by orthogonal polynomials (as q- 
Racah polynomials sit at the top of the q-Askey scheme - see |KS| ). Up to gauge equivalence, we obtain 
the weight of a horizontal lozenge with top corner (i, j) as: 

w(i, j) = Kg-' - (9) 

This weight was studied in |BGR10| . If we take k to or oo, we see the g-Racah weight is an 
interpolation between two types of weights: 

w{i,j) — q-' and w(i,j) ~ q^-' . 
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A direct alternative limit from the elliptic level is given by vi = V2 = p^f^,p — ?> (and then replace 
by (J or \/q). These two weights give rise to tilings weighted proportional to gVoiumc ^-Volume (^-^^j^gj-g 
Volume = number of 1 x 1 x 1 cubes in the stepped surface representing a tiling). This is the g-Hahn 
weight (as the analysis leads to g-Hahn orthogonal polynomials). The frozen boundary phenomenon for 
this type of weight was first studied in |KO07) , and then via alternative methods in [BGRlOj . 

Finally, the Racah weight is the limit g — >■ 1 in ([9| (we denote k = \ogq{K) and need k — )■ 1 as q — t- 1). 
The weight function becomes 

w{i,j) = k + j. 

Notice in all these limits the weight of a horizontal lozenge is independent of the horizontal coordinate 
of its top vertex. 

Taking these limits corresponds to the hypergeometric hierarchy of special functions involved in the 
analysis via the orthogonal polynomial (OP) or biorthogonal elliptic functions (down arrows denote 
limits) : 

Elliptic hypergeometric (elliptic weights; elliptic biorthogonal ensembles) 
g-hypergeometric (g-weights; q-OP ensembles) 

a- 

Hypergeometric (uniform/Racah weight; Hahn/Racah OP ensemble) 

As a side final note, the most general degeneration of the weight is the top level trigonometric limit 
p ~> 0, which gives rise to a 3 parameter family of weights (the use of the word trigonometric here should 



not be confused with its usage in Section 2.3). Being more general (more parameters) than the g- Racah 
limit, its analysis requires q rational biorthogonal functions rather than orthogonal polynomials. We will 
not use this limit hereinafter, as we can approximate the trigonometric level by choosing p really small 
at the elliptic level. 



2.5 Canonical coordinates 

It will be convenient for various computations to express the geometry of an elliptic lozenge tiling 
in terms of coordinates on a certain product of elliptic curves. First we will introduce 6 parameters 
A, B, C, D, E, F depending on q,t, S,T, N,vi,V2 (note we have listed, other than q, 6 parameters, of 
which 4 are discrete and dictate the geometry: t,S,T,N). t here is a (discrete) time parameter and 
ranges from to T. It will be explained better in Section |3] It corresponds to the fact that we will be 
interested in distributions of particles on a certain vertical line: that is, tilings of hexagons that have 
prescribed positions of particles (or holes) on the vertical line with horizontal coordinate t. The set of 
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parameters is: 



B = q*/2+S/2+T+l/2 



^ t/2-S/2-7V+l/2 L 



V"l^^2 
^ ^-t/2+S/2-Af+l/2 1 



(10) 



V^lW2 



^ ^ ^-t/2-s/2+i/2 Z:^ 



Observe that 



q^^-^ABCDEF = q. (11) 



Recall that the weight function (to be more precise, the ratio of weights of a full to an empty 1x1x1 
box - see ([s])) depends on the geometry of the hexagon via the three parameters ui,U2,U3 (J^Ufe — 1) 
which in the coordinates are: 



Us = q^^^^^ /viV2 

What we want is to change coordinates from (?,j) (2 dimensional) or (x,y,z) (3 dimensional) to 
(ui, 'U2, wa) via the above formula. We call these new coordinates canonical. In practice each line of 
interest in the geometry has an equation in the plane which can then be translated in terms of 



the Mfe's by solving in ( 10 ) for t, S, T, N, vi,V2 in terms of A, B, C, D, E, F. We thus find the following 



equations for the relevant edges of our hexagon: 
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1 /2 

Left vertical edge (corresp. equation : i — 0) : ^ — q~^vi/v2 — ( ^ ) E^q~^^^ 

1 /2 

Right vertical edge (corresp. equation : i — Tj : — — q vi/v2 — B q ' 

U2 \ DBF ) 



NW edge (corresp. equation :j = «/2 + iV): — = (7 



SE edge (corresp. equation : j = i/2 - - S*)) : ^ = q^^""^ 

Ml 



/^SC\^/^^_3 (12) 



f ABCV'"^ 



NE edge (corresp. equation : j ^ -i/2 + S + N) -.^ ^ q'^^+^^vivj = 

-"3 

SW edge (corresp. equation : j = —il2) : — = q^^viv^ = \ ) J:'"q 

U3 \DBFJ 

Vertical particle line (corresp. equation : i — t) : — — q~^*~^ 

U2 



V1/V2 



\DBF J 
DBF 



q 



ABC 



Remark 2.4. We can see from above that there exists a bijection between the six bounding edges of 
our hexagon and the 6 parameters A,B,C,D,B,F. That is, to an edge we assign the parameter that 
appears to the power ±3 above. The 6 parameters are not independent (they satisfy one balancing 
condition ABCDBF = q'^~'^^), but then neither are the 6 edges (they must satisfy the condition that 
the hexagon they form is tillable by the three types of rhombi, which in this case tautologically means 
the edges form the 6 visible frame-edges of a rectangular parallelepiped) . See Figure [7j 




E 



B 




Figure 7: Correspondence between edges and the 6 parameters. 



With ( 12 ) in mind we have a (local) map — ?> (where is isomorphic to the subvariety of E^ 
with coordinates (ui,M2,M3) and relation Yi^i — 1) which embeds our hexagon in E^: 



ihj) iui,U2,U3)- 
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Note that is the square of a real elhptic curve if parameters are chosen so that the weight ratio 
(of full to empty 1x1x1 box) is real positive. Hence as E is homeomorphic to a circle or a disjoint 
union of two circles, the above embeds our hexagon in a 2-dimensional real torus (base field = K). 



3 Distributions and transition probabilities 

In this section we compute the A'^-point correlation function and transitional probabilities for the model 
under study. We refer the reader to the Appendix of |BGR10j for the relevant application of Kasteleyn's 
theorem which makes these computations easy and to Kasteleyn's original paper for the theory itself 
( Kas6X). 

Take a collection of N non-intersecting lattice paths in n{N,S,T). Fix a vertical line inside the 
hexagon with horizontal integer coordinate t (0 < t < T). This vertical line will contain TV particles 
X — {xi < ... < xn) G J^v't- Depending on the geometry of our hexagon, there are four ways in which 
we can fix a vertical line with horizontal coordinate t inside a collection of N non-intersecting paths in 
Q{N, S,T). They are described below (see also Figure [s] in which the four cases are depicted - we only 
depict the outside bounding hexagon and the middle vertical line which is the desired particle line): 



Case l.t<S, t<T-S, 0<Xk<t + N ~1 

Case 2. S <t <T - S, 0<Xk<S + N-l 

Case 3. T-5<t<S', t + S - T < Xk < t + N - 1 

Case A.t>T-S, t>S, t + S - T < Xk < S + N - 1 




Figure 8: The four ways of choosing a vertical particle line (dotted) inside a hexagon. In all cases N — 5 
particles, T — 8, S €z {3, 5}. The middle vertical line in any hexagon is the particle line. 

We make the following notations: 

Lt(X) = sum of products of weights corresponding to holes (horizontal lozenges) to the left of the 
vertical line with coordinate t. The sum is taken over all possible ways of tiling the region to the 
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left of this line. Equivalently, it is taken over all families of paths starting at ((0, 0), (0, N — 1)) 
and ending at {(t, xi), {t, xn)). 



Rt{X) = sum of products of weights corresponding to holes to the right of the vertical line with 
coordinate t. The sum is taken over all possible ways of tiling the region to the right of this line. 
Equivalently, it is taken over all families of paths starting at {{t,xi), ...,(t,XN)) and ending at 
((T,5),...,(r,5 + iV-l)). 

Ct{X) = product of weights corresponding to the holes on this vertical line. 

Let 

'Pt,s{xk,Xi) = q--'9piq-'--',q-''+-'+^-'-^ViV2). (14) 

Remark 3.1. As observed in the introduction, Lpt,s{x, y) = —ipt,siyj x) so the product Y[k<i Vt,sixk, xi) 
is the "elliptic" analogue of the Vandermonde product (^*; ~ xi) (to which it tends in the limit 

p — > 0, g — 1 as explained in the Introduction). 

Proposition 3.2. We have 



Lt{X = {xi, ...,xn)) = const ■ ipt^sixk, xi)x 



k<l 



n 



Proof. This follows from an elaborate calculation and Lemma 10.2 in Appendix A of [BGRIO] (which 
is in essence an application of Kasteleyn's theorem and the computation of the appropriate inverse 
Kasteleyn matrix). 

First, as is in the case of the aforementioned lemma, we restrict ourselves to the case S < t < T — S 
(Case 2. in ( [l3| ; Computations are similar for the other 3 cases). Note in such a case we have N particles 
and 5* holes on the line with abscissa t. We then apply the particle-hole involution (as the weight in 
Lemma 10.2 in Appendix A of |BGR10| is given in terms of the positions of the holes ~ horizontal 
lozenges on the i-line). There are two types of products appearing in the total weight in question: a 
univariate one over the holes and a bivariate Vandermonde- like (again over the holes). For the first 
product, we just reciprocate to turn it into a product over particles (as the total product over holes 
and particles of the functions involved is a constant dependent only on t, S,T, N,q,p,vi,V2)- For the 
Vandermonde-like product, we note for a function / satisfying f{yi,yj) = —fiyjjyi) we have (up to a 
possible sign not depending on holes or particles): 

n fiyi'yj)= n f{xi,xj)x Yl I{u,v)y. 

l<i<j<S l<i<j<N 0<u<v<S+N-l 

l<i<N no<M<2;i fiXi,u) Y\xi<u<S+N-l /("' ^i) 



n 



where j/'s represent locations of holes (top vertices of horizontal lozenges) and x's locations of particles. 



We take / = (pt^s as defined in (14). Finally, in Appendix A of |BGR10| . the convention is that particles 
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and holes are counted from the top going down. This is opposite to the convention in this paper, so 
we substitute xi t-^ S + N — 1 ~ xi. After standard manipulations with theta-Pochhammer symbols we 
arrive at the desired result. □ 

Proposition 3.3. We have 

Rt{X = {xi,...,xn)) = const ■ ipt,s{xk,xi)x 

k<l 

n (n^-i^-s „-2t-S„, „1+T„, „1-T„, . „\ 

Proof. Similar to the previous proof except we use Lemma 10.3 in Appendix A of |BGR10| . □ 
Proposition 3.4. We have 

Ct{X = [xi,...,xm)) = const- || 



l<k<N 



g^'6'p(g2^' + l-*-'S'i;i-y2) 



Proof. This weight is (up to a constant not depending on holes or particles) the reciprocal of the total 
weight of the S holes (horizontal lozenges) on the t line and the latter is readily computed from the 
definition (|6|. □ 

Theorem 3.5. 

Proh{X{t) - {x^,...,xn)) = const ■Wi^t^sixk.xi))^ x [] g(2Ar-i).,g^(^2.,+i-t-s^^^^) ^ 



n 



k<l l<k<N 
^^^^^ dp{q, q^-S-t+T, q^-t-T-Sy^^qy^^ql+N-Sy^^^^ q^+^-ty^y^- q)^^ 

const -xiM^k,..)?- n q''^~'^^^oM^-F^) " 

k<i i<k<N ^pyq^q F^^i F^'i F^^'^k 



Proof. 

Prob{X{t) = {x^,...,XM))'^Lt{X)Ct{X)Rt{X). 



□ 



Remark 3.6. The above distribution is what was called in the Introduction the discrete elliptic Selberg 
density. That is to say, 

Prob{X{t) = {x^,...,Xn)) = ^x{q^^-^F^\q'^ ,q''-^AF,q^-\pB)F,q^~^CF,q''-^DF,q''-^EF) (17) 

where A e m" (to — S + N — l,n = N) and Xi + N — i — xjv+i-i (to account for the fact that 
Xi < X2 < ■■■ < xpf whereas partitions are always listed in non-increasing order). The particle-hole 



involution invoked in Proposition 3.2 then takes the following form. If Xp is the partition associated to 



19 



the particle positions (at time t) via the above equation and Xh is the partition associated to the whole 
positions at the same time (in the case above, there are S holes), then 

Ah = (Ap)' 

where A'^ denotes the complemented partition corresponding to A e m" (A^ — m — A„+i_i) and A' 
denotes the dual (transposed) partition (A^ = number of parts of A that are > i). The fact that both 
probabilities (in terms of holes and in terms of particles) are A-symbols can be observed directly as 
shown in Proposition 3.2 or using the following relations appearing in |Rai06] : 

Ay{a\...h-;l/q) = Axia/q^\...b,...;q) 



We will for brevity denote the measure described in Theorem (3.5 1 by pg t (note it also depends on 
N, T, vi, V2,p, q, but it is the dependence on S and t that will be of most interest to us). Observe we can 
transform the factor 

appearing in the univariate product of the above probability into something proportional to 

9p(g^-*-V,g^+^t;2) 1 



(g2-JV-t-S-T^^^ ^2-JV^2) ep{q^+'^-t-Svi,q-^+t+S+T/vi,q^+^+Tv2,q-'=/v2)N-i 



by using 



^JV-l.„^ _0p{A;q),9p{Aq-;q)N-i , , _ q^^'^'^-epiA; q),9p{q/A; q) 



0p[A;q)N^i ep{q^ ^/A;g)Ar_i 

and absorbing into the initial constant anything independent of x (of the particle positions Xk)- After 
using ( 10 ) our probability distribution becomes 

Proh(X{t) = {xi,...,xn)) = 

const ■Xliv.A-^^^x,)?^ n ep{B{Fq^^)^\l{Fq^^)^^-q)^_, ^ 11 ^'^^ 

where 



q-0piq^-+'-'-^'v,V2)9p{q'-^-\q'-^-S, q^-'-^ v,, q^+^V2, q'~^ViV2, q'-'-'viV2; g) 
9p{q^''~'^ViV2)ep{q, q^-S-t+T^ q^-^-t-T-S^^^ q2-N ql+N-S^^^^^ gl+^-*l^l«2 ; g). 

q-epiF^q^-)0piAF,BF{^j^^)\CF,DREF{^^^^^ 

0,{F')Op{^q, h (^^M) ' , gg, g<7, f g (^iM:) ^ , q; q)^ 
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w is the weight function for the discrete eUiptic univariate biorthogonal functions discovered by 
Spiridonov and Zhedanov (see [SZOOj . |SZ01) ). It is of course also the discrete elhptic Selberg density for 
= 1 (hence a A-symbol in n = 1 variable as seen in Q). Notice in ( 18 1 above B and E play a special 
role, as does F. This will become more transparent in Section [g] w is elliptic in q, ui,W2 and qi^^^^"^^^} 
(or, analogously, in A, B, C, D, E, F, q). 

Remark 3.7. Note that in the definition of w above, the first line is given in terms of the geometry 
of the hexagon and the choice of the particular particle fine (Case 2. in (13) as previously discussed), 
while the second line is intrinsic and the geometry of the hexagon only comes in after using (10). We 
can also define the equivalent of ( 10 ) in the other 3 cases described in ( 13 ) (and the three other choices 
of 6 parameters differ from (10) by (a): interchanging S ant t, (b): shifting the 6 parameters in (10) 
by or (c): a combination of both (a) and (b)). We will not use this any further, as all 

calculations will be done in Case 2. from (13). 

Remark 3.8. The limit vi = V2 = Ky^, p ^ gives the distributions present in |BGR10| at the g-Racah 
level. Also, as will be seen in Section[6j such probabilities are structurally a product of a "Vandermonde- 
like" determinant squared (the first two products in (18)) and a product over the particles of univariate 
weights of elliptic biorthogonal functions. Indeed, under the appropriate limits, one can arrive from ( 18 ) 
to a much simpler (prototypical) such iV-point function: the joint density of the eigenvalues of a QUE 
N X N random matrix. 

The transition and co-transition probabilities for the Markov chain X(t) are given by the next two 
statements. 

Theorem 3.9. If Y = {yi, i/n) and X = (xi, x^) such that yk — Xk £ {0, 1} Vfc, then 

ft+i,s{yk,yi) 



Proh{X{t + 1) = Y\X{t) =X)= const ■ fT x TT w.ix^) TT ^«o(^fe) 

'■PtS\Xk,Xi) 
k<l ' k:yk=Xk+l k:yk=Xk 



where 



g""gp(g"+i^^-^, g"+^+it;2, g""^+^i'it;2) 

""'^^^ " 0p(g2-+i-*-5i;i«2) 

Proof. The formula 

^ ^ ' ^ ^' ' Lt{X)Ct{X)Rt{X) 

^ Ct+i{Y)Rt+^{Y) 
Rt{X) 

along with the formulas for L, R and C yield the result. □ 
Theorem 3.10. If Y = {yi, y^r) and X — (xi, xn) such that yk — Xk € {0, —1} Vfc, then 

Prob{X{t - 1) = Y\X{t) =X)= const ■ TT "Pi-^^S^y^^Vi) ^ TT ^(^^.) TT ^^(^^) 

k<l • k-.yk^Xk-l k:yk=Xk 
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whe 



w'oix) 
w[{x) 



x-N-t+l 



-t-S+l 



V1V2 



V1V2) 



Proof. 



Prob{X{t - 1) = Y\X{t) =X) = 



Lt_,iX)Ct-iiX)Ct{Y)RtiY) 
Lt{X)Ct{X)Rt{X) 

Lt-i{Y)Ct-i(Y) 
Lt{X) 



□ 



We are now in a position to define six stochastic matrices (Markov chains) needed in what will follow. 
Their stochasticity along with other properties will be proven in Section |4j although we know the first 
two are stochastic as they represent the transition probabilities obtained in this section. To condense 
notation, we denote 



Let: 



t+ 



s+ 

:,S.t 



t-P. 



s- 



[0,1] 
[0,1] 
[0,1] 
[0,1] 
[0,1] 
[0,1] 



be defined by: 



const ■ u.^, "P'+'fy^^y') . n 



(yk,yi) rr g-"fcep(Azfc,Bzfc,Czfc,9i-"zfc/ABC) 



k<l i^t_s(a;fc,x,) llk:yk=Xk + l 6p{z? 



X 



Pt^iX,Y)=^ n,,^^,,^ .-^-"^^^(..M.W|,^./c.,"-^..ABC) (19) 
0, otherwise 



Vt-i, siVk.yi) T-f g~"'°~"+'gp(zfc/-P,Zfc/-E,Zfc/f.g"~'zfcg-EF) 



const ■ n.^, '^'-'■f^y'-y'' ■ n 



X 



0, otherwise 
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const ■ U^^, vt.s+ivv,,yi, . rj 



X 



0, otherwise 



Vt,s-i{yk,yi) rr g-="'=ep(B2:fc,Czfc,F^fc,gi-"2:fc/BCP') , 



const ■ U,^, vt,s-iyy,,yu . -rr . x 

0, otherwise 



n yt,s+ife.i/') n q-^k-«+Hp{zk/D.Zk/E.Zk/A.,q''-^ZkDEA) 

const -lY^^i ^^ .,(^^^^,) • llfe:y,=,;,-i e^(^2) X 

0, otherwise 

(23) 



COTISI j^j.^; v't,s(a:fe,a:i) \.\.k:Vk=Xk-l 0p(zl) 

0, otherwise 

. . .(^^^ 

The normahzing constants are independent of the x^s and the j/fc's. They wih become exphcit in 
Section IH 

Note that t-Pg-i under interchanging t and 5", becomes Pt-^ ■ Under the same procedure t+Ps+ 
becomes Pf_^ . We can think of PtX (Pt^S) as a Markov chain that increases (decreases) t, while t±Ps+ 
{t±Ps-) increases (decreases) S. 

Remark 3.11. In the q-Racah hmit vi = V2 = K^/p, p 0, the chains t±Ps+ coalesce into one (-Pj^ 

s,t 



in jB(;Rin) ;i. Likewise for t±P|l* 



4 Elliptic difference operators 

In this section we explain how recent results on elliptic special functions and elliptic difference operators 
intrinsically capture the model we described thus far. The main two references are |RailO| and |Rai06| 
and we will state results from these without going into the proofs (with a few exceptions where the 
proofs are short and revealing of common techniques employed in the area). The focus will be on 
certain elliptic difference operators satisfying normalization, quasi-commutation and quasi-adjointness 
relations. We define them abstractly in the first subsection. We then turn to motivating the definitions 
and interpreting the operators probabilistically. 
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4.1 Definitions and some properties 



In [RailO] (see also |Rai06j for an algebraic description) Rains has introduced a family of difference 
operators acting nicely on various classes of i3C„-symmetric functions. To define it, we let tq, ri, r2, r^, £ 
C* satisfy rorir2r3 —pq^^". Then define ©(r'o, J"!, ?'2, ''s) (also depending on q,p,n) by: 



(2)(ro,ri,r2,r3)/)(...Zfe...) = 2^ [[ " 2a,. [[ -IfTZi^TZ^fi-^ ^ ^k-)- {25) 

<Te{±l}" l<fc<n '^P^^k ) l<k<l<n k ) 

Remark 4.1. The difference operator above described is the special case t — q oi the more general 
elliptic {q, t) difference operator mentioned in the references. 

In view ofrQrir2r^ — pq^^" we will break symmetry and denote the difference operator by H{rQ, ri,r2), 
the fourth parameter being implicit from the balancing condition. 

Remark 4.2. © takes i?C,i-symmctric functions to i3C„-symmetric functions. 

By letting ® act on the function / = 1, we obtain the following important lemma, whose proof we 
sketch following |RailO| : 

Lemma 4.3. For rorir2r3 — pq^~" we have 

ET-r UoKsO^pi'^s^k'') T-r f^pi.l^k'' ^V) Tl a ( k k k \ 

11 2..^ 11 y7^^V^= 11 Opil r^r^^<l ror2,q r,r2) 

o-e{±l}" l<fc<ri "py'^k > l<k<l<n P^ k II 0<k<n 

Proof. By direct computation the LHS above is invariant under Zk — pzk for all k (this is insured by 
the fact rQrir2rj, — pq^^"). It is also _BC„-symmetric (invariant under permutations of Zk^s and under 
Zk — >■ Finally, by multiplying LHS by i? = Yik ^k^^pi'^k)Y[k<i "^(^k, zi) we will have cleared 

potential poles of the LHS. Because R is i3C„-antisymmetric the result will end up being a multiple of 
R: R ■ LHS = const ■ R showing LHS has no singularities in the variables and is thus independent of the 
Zi's. Evaluating then at Zi — rgg""' yields the result. Observe the main point here was to prove the 
LHS is elliptic and has no poles in the variables, and indeed any analysis that shows this will prove the 
result. □ 

Hereinafter we will use 2) for the "normalized" difference operator (so that ©(rp, ri, r2)l = 1) follow- 
ing Lemma |4.3| 

The difference operators described above satisfy a number of identities, including a series of com- 
mutation relations. For an elegant proof which relies on the action of these operators on a suitably 
large space of functions (more precisely, on the action of the difference operators on _BC„-symmetric 
interpolation abelian functions), see [RailO] or |Rai06| . 

Lemma 4.4. // U, V, W, Z are 4 parameters, then 

'D{U, V, W)D{q^/^U, q^/^V, q^^^^Z) = ©(t/, V, Z)D{q^l'^U, q^''^V, q-^''^W) 

Next we look at the action of the difference operators on special classes of functions. For A e m" a 
partition, let 
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, , T-r ni<;<r»+n^p(V '^fc ) 

i<fe<„ lli<;<n "^pl^y -^fc 7 

By direct computation, we see that dx{...uq^''^"^'' ...) = S\^^c\. 

Remark 4.5. d\ is a special version of the interpolation theta functions P^^"^'"'\...Xk-.-; a, b; q;p) defined 
in |Rai06| (matching the notation in the reference with ours, a = u,b = g~"'~"+^/a)). They are defined 
(up to normalization) by two properties: being i?C„-symmetric of degree m (which happens for d\a) 
and vanishing at fi ^ X (which trivially happens in our case). 

If we now define \ = — we see that 

c)A(...ug^'=+"-^..) = ^A,M (26) 

so in a precise way, dx is an interpolation Kronecker-delta theta-function. We then immediately have 
the following proposition: 

Proposition 4.6. Fix t e {±1}". Let Zk = uq^>'+"-''. Then 



Proof. Immediate by substituting into the definition of the difference operator (25). For any a ^ t, 
qak/2-Tk/2 will be of the form with /i 7^ A and the corresponding summand will be 0. □ 

A useful final property of the difference operators is their quasi-adjointness. It was shown in [RailO] 
that the ©'s satisfy a certain "adjointness" relation that wc will need in the next section. We start with 
6 parameters io, ii, ^2, ^3, uo, ui satisfying the balancing condition 

q'^"'^^tQtit2t3UoUi = pq. 

We fix the number of variables at n and A will be a partition in to" . As in the introduction, we 
denote li = Xi + n — i. We define the discrete Selberg inner product (, ) (depending on p, q and the 6 
parameters) by 

(/'5> - I E /(■•■^o9'^.05(■••^o9'^.0AA(g'"-'^^|g^g"-'^otl,g""'^o^2,g"-Ho^3,g""'^o"o,g""'^oWl;<z) 

ACm" 

(27) 

where /, g belong to some sufficiently nice set of functions (we will assume they are _BC„-symmetric) 
and Z is an explicit constant that makes (1, 1) = 1. This is a discrete analogue of the continuous inner 
product introduced in jRailO] and can be obtained from that by residue calculus. 
If the above conditions are satisfied, then f |RailOj ): 

(©(wo, to, h)f, g) = (/, ©(u'l, 4, t'^)g)' (28) 

where 

(i'o, AAA. < O = ('Z'^'io, gi/'ti, g-^/'t2, 9"'/'i3, 9'/'wo, q-^'''ui) 
and (, )' is the inner product defined in \27\ with primed parameters inserted throughout. 
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4.2 Interpretation of difference operators and their properties 



We now show how the difference operators and their properties discussed in the previous section can be 
given probabihstic interpretations. 



Remark 4.7. Observe from (10) that q^'^-^ABCDEF = 1 



In what follows hk {h'j^) is the location of the fc-th particle on the vertical line i = t [i = t -\- 1) m 
the (i, j) frame (note according to the f — )■ t + 1 dynamics the particles move either up or down by 1/2). 
We can prove the following proposition: 



Proposition 4.8. For A,B,C,D,E,F and Zk — Fq given by (10), the summands in 



(©(AB,C)l)(...Zfe...) 



(see (|25[) ), appropriately normalized using (4.3), are equal to the transition probabilities (entries in the 



stochastic matrix) P^_^_ {H,H') defined in (19) (after switching coordinates from {x,y) back to {i,j)). 
This statement also holds for: 



'D{D,E,F) and P(_: 
'D{A,B,D) and t+P. 



'D{B,C,F) and t+Pg 
'D{D,E,A) and t^P. 



S,t 
S,t 



S,t 

s+ 

S,t 



'D{E,F,C) and t^Pg 

Proof. 1 will only prove the statement for 'D{A, B,C) and t+ (the equivalent statement for 'D{D, E, F) 
and t— is proved much the same way). The proof is immediate in view of (10), the change of variables 



{X,Y) M- {H,H') in (19) (to the (i,j) coordinates) and the following observations. 



First, a choice of e {±1} for all k in the definition of 'D{A, B,C) is equivalent to a choice of 
which particles move up/down from the position vector H (at vertical line t) to the position vector 
H' (at vertical line t + 1). If — 1, the corresponding fc-th particle at vertical position h^. moves 
up to h'f, = hk + 1/2 (and if dj. = —1, the fc-th particle moves down). Next observe that in the 
univariate product appearing in any term of {(D^A, B, C)l)(...Zfc...), we can change 9p[uz~^) (6 = 1, 2) to 
9p{z\/u) by the reflection formula for theta functions and it will now match with the univariate product 
appearing in Pi^ . The product Y\k-y^=Xk+i^"-^^yk=xS'"^ '^'^^ indeed is identical (modulo constants 
independent of the particle positions) to Ofc h' =hk+i/2i---) Y\k h' =/ifc-i/2(---) which is nothing more than 

n.:.,=i(-m.:.,=-i(-)in(iii|. \ . 

The elliptic Vandermonde product Jlfc<i appearing in ( 19 ) is the same product (modulo constants 
independent of the particles) as the Vandermonde-like product in any term of {'D{A,B,C)l){...Zk-..) 
once we've transformed (in the latter product) Op{zi/zk) into dp{zk/zi) and 9p{l/zkZi) into 0p{zkZi) 
(picking up appropriate multipliers in front that will be powers of q appearing the Vandermonde-like 



product in (19)). The extra powers of q appearing in (19) will also surface in the difference operator 
once we've performed the aforementioned transformations. Finally observe that the ratio '^*+i'S(^fc:fei) 



Vt,s{hk,hi) 



reduces (modulo the power of q up front already accounted for) to a ratio of only 2 theta functions (of 
the 4 initially present) because either h'^. — h\ — hk~hi or 



-h[ = hk- 



k and I moved both in the same or in different directions) 



- hi (depending whether particles 

□ 
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Remark 4.9. We describe how the difference operators capture the particle interpretation of the model 
intrinsically. In their definition specialized appropriately as in the statement of the above proposition, if 
two consecutive particles k,k + 1 are 1 unit apart (hk+i — hk = 1), the bottom one cannot move up and 
the top one down to collide because the summand in the difference operator is (indeed Op{qzi~z'j^^-^) = 
9p{l) = in the cross terms). Thus, the non-intersecting condition on the paths is intrinsically built 
into the difference operator. A similar reasoning shows that top-most and bottom-most particles are 
not allowed to leave the bounding hexagon either. To exemplify, for the difference operator 'D{A, B, C) 
corresponding to the t ^ t + 1 transition (particles moving from left most vertical line to the right), we 
observe that the restriction on top (bottom) particle is not to cross the NE (SE) edge labeled C {A) in 
Figure [t] (or indeed not to "walk too far" to the right by crossing the B edge). However A and C are 
two of the parameters of the difference operator, and the corresponding terms in the univariate product 



in the appropriate summand in ( 25 ) become once the top (bottom) particle tries to leave the hexagon. 
Same reasoning applies to the particles not being able to "walk too far right". Hence the difference 
operators intrinsically capture the boundary constraints of our model. 

Remark 4.10. Proposition 4.8 is even more general, as we obtain (g) — 20 different stochastic matrices 



(Markov chains) from the 20 different difference operators (6 of them already described). 

We are now in a position to prove that the 6 matrices defined in section [3] are indeed stochastic and 
measure preserving. 



Theorem 4.11. 



Y 
Y 



Ps,t±iiY) ^J2p'±\X,Y) ■ psAX) 

X 

Ps±i,t{Y) = E t±Ps±iX, Y) ■ psAX) 



X 

Proof. There is one way to prove these statements which works for 4 of the 6 matrices. Observe that 
the results for t± follow from Theorems |3.9| and |3.10[ and then to observe that under t -fr^ S, we have 

and then under interchanging S and t, P^_^ becomes t+Ps+ P^Z becomes t-PsL, respectively). 

This idea worked both at the q-Racah level and Hahn level (see [BGRlOj and |BG09| ). 



Alternatively we ca n ob serve that the first two equalities are, by using ( 10 ) and Proposition |4.8[ 
restatements of Lemma 



{D,E,F) (for Pt'), {A,B,D) (for t+P|;'), (B,C,F) (for *+P|l^), {D,E,A) (for t_P|;^), (P, F, C) (for 
4_P^1*). Moreover, the normalizing constants that we omitted in defining the transition matrices can be 
recovered easily from Proposition |4.8| 

The last two statements are a special case of the adjointness relation. We will prove the third 
statement for the t+ operator. Similar results exist for the other 5 operators. We recall that ps,t{X) is 
nothing more than the discrete elliptic Selberg density 

A,, (g^^-^F^Ig^, q^'-'AF, q''-\pB)F, q^'-'CF, q^'-'DF, q^-^EF) 



4.3 



for difference operators corresponding to parameters {A,B,C) (for P^^ 

\S,t\ ( T-) rp\ /£ T-)S,t\ / J-\ 771 A\ ( £ TjS ,t\ 
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defined in the introduction, with \x.k + n — k — Xn+i-k- We also define the partition Ay to be the one 
corresponding to vertical line t -\- 1 and particle positions given by Y: Xy^k + n — k — yn+i^k- Then one 
sees ps,t+i{Y) = Y.x Pt+i^^Y) ■ psA^) is equivalent to: 

(2)(AS,C)c)A,,l) = (0 A,, (29) 



where the prime parameters and (, )' are defined in the previous section. The above equality (29) is only 
"morally correct" as we encounter the following issue: the (summands in the) difference operators D 
correspond to transitional probabilities in the coordinates where particles move up or down by 1/2 



from the t vertical line to the t + 1 vertical line (from Proposition 4.8). Dx, Ax as well as the definition 



of the inner product (27) correspond to coordinates {x, t) where particles cither move horizontally 1 step 
to the right or diagonally up by 1 from vertical line t to vertical line t + 1 (see the previous subsection 
and recall +n— fc = Xn+i-k)- But this can be easily fixed since = {t,x — t/2). However, writing 
an "approximate" version conveys the meaning of the quasi-adjointness of the difference operators in a 
notationally uncluttered way. 

With the previous comment in mind, the right hand side in (29) equals Oa^, (•.•-Fq^''^" ''■•O^^t — 
A';^^ = ps^t+i{Y) (observe A' = A with prime parameters corresponds to the distribution of particles 

at the line t + 1) while the left hand side equals Y.Xx Prob{\Y\\x) ■ Aa^ = Ex Pt+{Y\X) ■ ps,t{X). 
The result follows. □ 

We finish this section with a graphical description of the 6 Markov processes described thus far. The 
key is to look at the domain and codomain of the difference operators in canonical coordinates. We will 
exemplify with the difference operator 'D{A, B, D), corresponding to Markov chain t+Ps+- Recall this 
Markov chain quasi-commutes with the t — > i + 1 chain. The key is the following relation (a restatement 



of Theorem 4.11 1 



Prob{Y\X; A, B, D)Prob{X; A, B, C, D, E, F) ^ Prob{Y; A', B' , C , D' , E' , F') where 

X 

{A',B',C',D',E',F') = {q^A,q^B,q-^C,q^^D,q-^E,q-^F) 



We note t+Ps+ corresponding to difference operator 'D{A, B, D) maps marked random tilings of 
hexagons determined by parameters {A, B,C, D, E, F) to random tilings of hexagons determined by 
parameters {A' , B', C , D' , E' , F') (marked here refers to the particle line corresponding to parameter t). 
We figure what happens to the edges of s uch hexagons when parameters get shifted by g^*^^/^ by using 



(12) (canonical coordinates). Figure 4.2 is a graphical description. In particular, we observe t+Ps+ 
increases 5 by 1. Similarly for the other difference operators: they increase (decrease) 5* or t by 1 while 
leaving the other constant. 



5 Perfect Markov chain sampling algorithm 

5.1 The S+1 step 

In this section, which follows closely the notation and proofs of IBG09) (see also [BGRlOp . we will define 
a stochastic matrix 

P^^s+i ■■ ^{N, S, T) X n{N, S+l,T)^ [0, 1] 
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Figure 9: Action of the difference operator 'D{A,B,D) on a tiling of a iV = 2, 5 = 4, T = 7 hexagon 
drawn in canonical coordinates. The source is marked 1 and the destination 2. Only edges relevant 
to the model are considered: the 6 bordering edges and the particle line at horizontal displacement t 
from the leftmost vertical edge. Note the slight shifting, the increase in S" by 1, and the fact that the 
particle line's displacement from the left vertical edge (= t) is kept constant (though particle positions 
are shifted by a third step). 



that is measure preserving: it preserves the elliptic measure /i(A^, 5, T) - the total mass of a hexagon 
tiling (collection of N non-intersecting lattice paths) in il{N, S,T). Recall /i is defined as a product of 
the weights of the individual horizontal lozenges inside the hexagon. Viewed as a Markov chain, the 
input for Ps^g^i is a hexagon of size a x b x c and the output a hexagon of size a x (6 — 1) x (c + 1). 
Both the input and the output will turn out to be distributed according to ii{N, S, T) and /i(-/V, 5 + 1, T) 
respectively. 

Given a collection of non-intersecting paths X = {X{Q), ...,X{T)) e rt{N, S,T), we will construct 
a (random) new collection Y = (F(0), F(T)) e il{N,S + l,r) by defining a stochastic transition 
matrix Pg^g_^_l{X, Y). Observe that F(0) € x^'^^'° = (0, - 1) is unambiguously defined. Next we 
perform the sequential (inductive) update. That is, the procedure which produces a random Y{t + 1) 
given knowledge of Y{0), Y{t) and X which we assume to have already been obtained. Y{t + 1) will 
be defined according to the distribution 



Pf+^^\Y(t), Z) ■ t+P^+^''+\Z, X(t + 1)) 
Prob{Yit + l) = Z)= ^ ,\ , - 

(pf/i'*.,+p|+i-*+i)(y(t),x(< + i)) 
^ ,_p|:,*+^(x(t + i),z) ■ pf_+'^'+\z,Yit)) 
(,_p|;*+i . Pl+'^'+'){X{t + l),Y{t)) 



(30) 



where the last equality follows from the fact that ps+i.t+i{A)Pt-^^'*^^ i^j B) = ps+i^t{B)Pfj^^'*{B,A) 
(this is nothing more than the equality Prob{An B) = Prob{A)Prob{B\A) = Prob{B)Prob{A\B)) . 
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We define the matrix Ps^s+i ■ ^{N, S, T) x n{N, S+l,T)^ [0, 1] by 



T-i pff'-'(y(t),y(t+i))-t+p|+'-'+'(y(t+i),x(f+i)) 

t=0 TTJH+TTt— ^^T+TTt+TT 



^s^s^i = S if UjJiPf-,^'' ■ t^P^^'-'^'){Y{t),X{t+l)) > (31) 
0, otherwise 

Theorem 5.1. The matrix Ps^s+i stochastic and ^-measure preserving, in the sense that 

^^{N,S+l,T)iY)= Ps^s+l{X,Y)^liN,S,T){X). (32) 

Proof, (following |BG09j ) We want to show that 

. ^ _ "tt' Ptl^'''iY(.t), Yjt + 1)) ■ ,+P^+'^'+\Yit + l),Xit + 1)) _ ^ 
^P.^.,,(X,y) = ^ [[ (P-M.^^p|+M.i)(^(,),^(, + ,)) - 

where the sum is taken over all Y = (Y{0), Y{T)) e n{N, S +1,T) such that 

T-l 

n {PtT'' ■ t+P^^'''^'){Y{t),X{t + 1)) > (33) 
t=o 

We first sum over Y{T) and because F(r) is distributed according to a singleton measure, the 
respective sum is 1. Next we deal with the sum 

pf/^-^-^(y(r - 2),Y{T - 1)) ■ ,+p|+^-^-^(y(r - 1), x(r - 1)) 

/pS+l,T-2 pS+l,T-l 



y(T-i) 



(P*/ • t+Ps- ' WiT - 2),XiT - 1)) 



over y(T- 1) satisfying {Pf+^'^~^ ■ t+P|+^'^)(r(r - 1), X(T)) > (because of (|33|)). Because of the 
quasi-commutation relations from Theorem |4.4| we have 

(pf++i'^-i • ,+p|+i'^)(y(T - i),x{T)) ^ (p|+i'^-i . ,+p|i^-i)(r(r - i),x(r)) 
> p|+^''^"^(r(r- i),x(p- i))P4'^'^"^(x(r- i),x{t)). 

We are summing over Y{T — 1) such that the LHS above is non-vanishing, but if it vanishes, then 
by the above inequality so does t+Ps-^''^~^{Y{T — l),X(r)) (one of the numerator terms in the sum 
over Y{T — 1) considered). This means we can drop the condition that (P^^ ' ^ • t+Ps~ ' ){Y{T — 
1),X{T)) > and sum over all Y{T — 1). We obtain 1 for this sum (the denominator is independent of 
the summation variable, and summing the numerator over Y{T — 1) we obtain the denominator). We 
next sum inductively over Y{T — 2) and so on until we are left over with a sum over Y{0). This sum 
only has 1 term, so we obtain the desired result. 

To show Psh^s+1 preserves the measure fi, observe first that 

M(iV,5,r)(X) = mo(X(0)) • Pf_fiX{0),X{l))...Pff-\XiT-l),XiT)) 
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where mg is the unique probability measure on any singleton set (in this case X '"). Then the RHS of 



( 32 1 becomes 



(34) 

Pulling out factors independent of the summation variables, replacing 1 = mo(X(0)) with 1 = mo(y(0)), 
using t+Pt''^{YiT),X{T)) = ,+P|+i^"(y(0), X(0)) = 1 and Pf+* • ,+P|l*+i = ,+P|l* • Pf^"^'*, we 
transform ( [34| into 



moiYio)) n p,r''(y(i),ni + 1)) X E n — . ^s,t o^+M 



Now we sum first over X(T), then over X{T — 1) and so on like in the previous argument to finally 
obtain on the LHS the desired result: 

T-l 

mo{Y{0)) n Pf+^''\Y{t),Y{t + 1)) = /i(Ar, S+1, T){Y). 
t=o 

□ 

5.2 Algorithmic description of the S" i— )■ + 1 step 

As before, whenever possible, we try to keep the notation similar to |BG09| . For a; G N we define 

9p{q^+\q^-^^-S-ivi,q^-S+T+iv2,q'=-^+^ViV2) ^ ep{q^^~-^-S-^ViV2y 

Note p also depends on 5, T, vi,V2, q,p, but we will omit these for simplicity of notation. Also note 
p is an elliptic function of q, q^ ,q^ ,q*,vi,V2,q^ ■ Consider (again omitting most parameter dependence) 



P{x;s) = Y[pix + i-l). 

1=1 

P is just a ratio of 5 length s theta-Pochhammer symbols over 5 others (multiplied by q^^^ to make 
everything elliptic). We define the following probability distribution on the set {0, 1, n}. 

Proh{s) = D{x- n){s) = .. . (35) 

For the exact sampling algorithm, given X — (X(0), X(T)) G il{N,S,T), we will construct 
Y = (y(0),...,r(T)) e n{N,S + l,T) by first observing that y(0) = (0, ...,A^- 1) is uniquely defined. 
We then perform T sequential updates. At step i + 1 we obtain Y{t + 1) based on Y{t) and X{t + 1). 
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Suppose X{t+ 1) = {xi,...,xn) e X'^'*+^ and Y{t) = (yi, yjv) G X'^+^'*. We want to define/sample 
Y{t+1) = (zi,...,ZAr) andX(t + l) satisfy - g {0,-1,1} (follows by construction 

from {Pf^^'* ■ t+Ps-^'^^^)(X{t), X{t + 1)) > 0). We thus have three cases, and we describe how to 
choose Zi in each: 

• Case 1. Consider all i such that Xi — yi = 1. Then Zi = Xi is forced. 

• Case 2. Consider all i such that Xi — yi = —1. Then Zi ~ is forced. 

• Case 3. For the remaining indices, group them in blocks and consider one such called a (fc, I) block 
(where k is the smallest particle location in the block, and I is the number of particles in the block) . 
That is, we have yi_i < fc — 1, yi+; > k + I and the block consists of 

Xi^yt^ k, Xj+i = yj+i = fc + 1, ...,Xi+i^i = yi+i-i = k + l -1. 

For each such block independently, we sample a random variable ^ according to the distribution 
D{k;l). We set Zi ~ Xi for the first ^ consecutive positions in the block, and we set Zi = Xi + l for 
the remainder of the / — ^ positions. We provide an example in Figure [l0| below: 

X(t) X(t+1) Y(t) Y(t+1) 

o o 
° o could not jump (first case) 



o 

- o — 



^ Split point determined by (S) 
T (block case) 

2 



o o 



s+1 



was forced to jump (second case) 



Figure 10: Sample block split. Same picture appears in |ljG09j for uniformly distributed tilings. 
Theorem 5.2. By constructing Y this way, we have simulated a S t-^ S + 1 step of the Markov chain 
Proof. We perform the following computation (and are interested in Case 3. described above, that is on 
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how to split a (fc, I) block; note Xi = Ui in the case of interest): 



P!+^'\Y{t), Z) ■ t+Pt'''+\z, Xjt + 1)) 

(pf/i'*.,+p™+^)(y(i),x(< + i)) 



Prob(Y{t + I) ~ Z) — — - — gqrj-^ s_|,i t+i-, ,^r, — ^ (factors independent of^) (36) 

-T-S-< 

W^-*-^«l«2) 

TT -,. gp(g^''^"^,g^'-^*-^-^^l,g^'+^+^^2,g^'-^+^^l^^2) .33. 



11 flr„2,.-t-s,„„ „A ^^^^ 



T-r gp(g-'-^-^, g-'-*-^-l^;l, g-'+^+iz;^, g"'~^-g-l^;l^2) 

11 ^^p(g2-.-*-s~l„^„2) ^ ^ 



a;. 



X n g---^^^^^^^^^ , , ^':ls^, ^ (40) 

We thus see the blocks split independently. The probability that the first j particles in a (fc, I) block 
stay put from Y{t) to Y{t + 1) (and the rest of Z — j jump by 1) is, by using the above formula: 



n 



i=0 

X 



^^(^2fe+2«-t-S-l^^„2) 



n 



Q^^q2k+2r-t-S+ly^y^^ 



X (factors independent of j) 



where in (36) we have gauged away everything independent of the split position j. This probability is 
nothing more than the distribution D we defined in ( 35 1 . This finishes the proof. □ 



5.3 Algorithmic description of the S ^ S — 1 step 

Similar to the Ps^s+i matrix described in the previous two sections, we can construct a Ps^s-i measure 
preserving Markov chain that takes random tilings in f2(iV, 5, T) and maps them to random tilings in 



f2(iV, S — 1, T). We proceed exactly as in Section 5.1 and will omit most details and theorems as they 



transfer verbatim from Section 5.1 and the previous section (we refer the reader to jBG09j for more 
details on this algorithm). Given X G fl{N,S,T) and Y{0),Y{1), ...,Y{t) already defined inductively, 
we choose Y{t + 1) from the distribution: 



We define 



P,l-'-\Y{t),Z) 



Prob{Y{t + 1) = Z) = 



P^I^^'+\Z,X{t + l)) 



S-l.t 

t+ 



S-l,t+l 



){Y{t),X{t + l)) 



(41) 



\Y{t),Y{t+l)) 



lt=0 



\Y{t+l),X(t+l)) 



PS^S-1=\ :f T-\T-lrpS-l,t S-l,t+l 

0, otherwise 



){Y{t),X{t + l))>0 



(42) 
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We will also sketch the algorithm for sampling using Pst-^s-i- We need to define the equivalent for 
p from the previous section. For x € N we define 

p [x) = — ^^^T^ rT, T^nrrr-, r X 



As before, p' is an elliptic in g, q^ , q^ , q* ,vi,V2,q^ ■ We also have P'{x; s) = 11;=! p'i^ + i — I) and 
the following distribution on {0, 1, ...,n}: 

Prob{s) = D'{x-n){s) = ^""'f .. . (43) 

Assuming we have X g ^1{N, S, T) with X{t + 1) = {xi < ... < xn) and inductively Y{0), ...,Y{t) = 
{yi < ... < t/Af), we sample Y{t + 1) (zi < ... < zjy) by first observing that Xi — yi€ {0, 1, 2} (because 
{Pfj^ ' ■ t+Ps+ ' )(^(0)^(^ + 1)) > 0) ^i^*^ then performing appropriate updates for the following 
three simple cases: 

• Case 1. For all i with Xi — ?/i = we set z,j = Xi. 

• Case 2. For all i with Xi — yi — 2 we set Zi = yi + 1. 

• Case 3. For the remaining indices (for which Xi ^ yi = 1), group them in blocks and consider one 
such called a (fc, I) block (where k is the smallest particle location in the block, and / is the number 
of particles in the block). That is, we have yi-i < k — 1, yi+i > k + I and the block consists of 



Xi =^ yi + 1 = k, Xj+i = yi+i + 1 = fc + 1, x,,+i-i = y.+i-i + 1 = fc + / - 1. 

For each such block independently, we sample a random variable ^ according to the distribution 
D'{k; I). We set Zi — yi for the first ^ consecutive positions in the block, and we set Zi = y^ + 1 for 
the remainder of the I — ^ positions. See Figure \T0\ 

An analogous of Theorem |5.2| exists and is proved in a similar way to show the above 3 steps are all 
that is necessary to simulate the Markov chain Ps^s-i- 



6 Correlation kernel and determinant al representations 

In this section we will show the process X{t) corresponding to a tiling of the hexagon is determinantal 
with correlation kernel given in terms of elliptic biorthogonal functions due to Spiridonov and Zhedanov 
f jSZ00| . [Rai06p . We start by a brief overview of the necessary facts about biorthogonal functions, and 
continue with the heart of the proof: an application of the Eynard-Mehta theorem. 

6.1 A brief overview of elliptic biorthogonal functions 

We will first gather together a few results about univariate discrete elliptic biorthogonal functions. The 
notation and exposition will mostly be following |Rai06j . We will need to make brief use of univariate 
interpolation abelian functions. They were introduced in [RailOj (see also |Rai06) for a description closer 
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to our purposes). They are, for a fixed integer /, BCi -symmetric (i.e., symmetric under x i— ^ 1/^) ratios 
of i?Ci-symmetric theta functions of degree I with prescribed poles and zeros. To wit: 



Observe i?;* has zeros at finitely many g-shifts of a and poles at finitely many g-shifts of h (up to 
taking reciprocals and shifting by p). The biorthogonal functions Ri{x;tQ : ^i, ^3; up, Ui) discovered 
by Spiridonov and Zhedanov f |SZOO| in the univariate case; see |RailOj for continuous and |Rai06) for 
discrete multivariate analogs) can be defined in terms of the interpolation functions as follows ( |Rai06) ) . 
Fix \p\ < 1, q as well as six parameters to, ^1,^2,^3, uq, ui such that totit2t3UoUi — pq. Then (dependence 
on p, q implied but not written) 

Ri{x;to : ti,t2,t3;uo,Ui) — (ifei?^(a;; to, Ug) = dji?;* + lower order terms 

o<fe<z 

where the formula for the 's is explicitly given in jRai06| and is independent of x (but of course depends 
on tQ,ti,t2,t3;uQ,ui,q,p and k). These functions have poles at shifts of Uq^ (we will say uq controls 
the poles of Ri{x;to : ii, ^3; "o, wi)). They are elliptic in the 6 parameters (provided the balancing 
condition is satisfied) as well as in the variable x. Furthermore, if in addition to the balancing condition, 
one also has 

toil = g-" (44) 

(for some m > an integer), the functions with poles controlled by uq and those with poles controlled 
by ui satisfy the following discrete biorthogonality relation on {0, m}: 

Ri{toq'';to : ti,t2,t3;uo,ui)Rk{toq'';to : ti,t2,t3;ui,uo)As{tl\q,toti,tot2,tot3,toUo,toUi) = Si^kQ 

0<s<m 

where As is the univariate weight discussed in the Introduction (also appearing in section |3| and 

ci = const ■ Ai{l/uoUi\q,toti,tot2,tot3, l/io^o, l/^owi)"^ = const ■ AiillqJoiiJohJoUJoUoJoUiy^. 

(45) 

The "hat" parameters are defined by the relations 

to - \/^^^, tot, = tot,, / = / (46) 

V pq to to 

for i = 1,2,3 and j = 0,1. The "hat" is an involution. Also observe the hat parameters satisfy the 
same balancing conditions the original parameters satisfy. They are important because by hatting we 
can exchange the variable and the index of the biorthogonal functions as follows: 

Ri{toq'';tQ : ti,t2,t3;uo,ui) = Rs{ioq^;io : ti, £2, 4; "o, Mi)- (47) 

The biorthogonal functions described above have to as a special normalization parameter (distin- 
guished among the t^'s). That is, Ri{to',to : ti, t2, ts; uq, mi) = 1. The normalized difference operators of 
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section |4] act on the biorthogonal functions as follows (note uo is special - it controls the poles, and to is 
also special as the choice of normalization): 

i)iuo,to,h)Ri{{q^/ha)q';q^/ho : q^'^h,q-^'H2,q-^'^h;q^'^Ua,q-^'^Ui) = Ri{toq';to : ii, ta, tg; Uq, Mi). 

(48) 

Finally, we can exchange to with another tj at the choice of picking up a factor (this is in essence a 
renornialization so that R takes value 1 at tj rather than tg): 

X Ri{x;to ■.ti,t2,t'i]Uo,ui) 
Ri(x;ti ■.to,t2,t3;uo,ui) = ——— — r (49) 

ni\tl',to ■ 11,12, Is; MOi'"lj 

6.2 Determinantal representations 

In this section we will show the processes t i-^ t ± 1 are determinantal point processes. For a review 
of such processes we direct the reader to jBorllj . We will do the calculation for the t i-^ t — 1 Markov 
process as it leads to less complicated formulas, but analogous results hold for i i-)- 1 + 1. 

For the remainder, it is now convenient to relabel and rescale the parameter set {A, B, C, D, E, F} 
as {^0,^1,^2,^3,^0,1*1} in order for certain symmetries to become more prominent (and in doing so, we 
will use the notation set forth in the previous section). To wit: 

A = t2, q^-^B = lii, C = t3, D = ti, q^-'^E = uo, F = to- (50) 

Note these parameters depend on t (the time parameter), and such dependence will be made more 
explicit when it becomes important. Notation is as in the previous section. Note uoMiio^i^2i3 — 1- Since 
the balancing condition for the biorthogonal functions requires a pq on the right hand side, we will again 
multiply Ml by p. These are the parameters of the univariate biorthogonal functions discussed in the 
previous section, uo and ui control the poles of the pair of biorthogonal functions. 

At the core of the computations will be the Eynard-Mehta theorem, which we now state in a 
"decreasing-time" form convenient for our computations (see |EM98) , |Borll| for a review and |BRQ5| 
for an elementary proof): 

Theorem 6.1. Assume we are given the following: 

• a discrete hiorthonormal system {fl,g\)i>o on 12{{0,1, L}) for each time t = 0, ...,T 

• a matrix 

Vt^,^^{x,y)=J2fr\tl\n9l(tlqy), 
l>0 

for n > 0, t — 1, T and a parameter to changing with time 

• a discrete time Markov chain X(t) (with time decreasing from T to 0) taking values in state spaces 
X* (set of possible particle positions at time t) with one dimensional distributions proportional to 

det (/fe_i(4'?'")) det {gLiitU"")) 

and transition probabilities proportional to 

deti<kj<N{vt^t-iixk,yi))deti<k,i<Niflz\(tl''^q'"')) 
deti<fe,,<Ar(/*_i(4(?-0) 
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Then 



Prob{xi £ X{ti), ...,Xs e X{ts)) = det {K{Tk,Xk;Ti,xi)) 

l<k,l<s 

where 



K{ti,Xi;T2,X2) 



-Y.s>Nf?{Q<f')9?{tl'<l^')^ if n < T2 



The first step in showing the required determinantal formulas needed to apply the Eynard-Mehta 
theorem is the following determinantal formula, a version of which was discovered by Warnaar (see 
[War02j Lemma 5.3 and Corollary 5.4 for comparison): 



Lemma 6.2. 



det i?;_i(zfc;io : <i, <2, <3; "cP^i) = const ■ TT(^(zfe,z/) TT 



where Zk — toq^'' , the constant is independent of the Zk 's and nonzero. 

Proof. This proof is essentially the same as that of Lemma 5.3 in |War02j . but is reproduced here for 
clarity. A first observation is that the constant in front of the right hand side will not matter much, 
and because it is ignored, the proof is somewhat simpler (of course, something has to be said about 
it not being 0). If we denote the left hand side by L and the right hand side by R, we notice both L 
and R are elliptic in the Zk^s (for R this is a direct calculation, and for L the biorthogonal functions 
inside the determinant are elliptic as mentioned in the previous section though one can just see this 
from the definition in terms of abelian interpolation functions). Fixing a variable Zk, we see poles for 
L/R come from the zeros of R or the poles of L. For the latter, the poles are controlled by uq but 
are exactly canceled by the zeros of l/i? appearing in the univariate product (one can see this from 
the definition of biorthogonal functions in terms of abelian interpolation functions) . For the former the 
zeros of R possibly leading to poles are Zk = zi, Zk = 1/ Zi ior I k (and p shifts thereof). Plugging in 
Zk = zi into L makes two columns the same, so L vanishes. Since univariate biorthogonal functions are 
-BCi-symmetric in the variable (a fact made explicit in the previous section in the definition in terms of 
abelian interpolation functions; see also |Rai06| ) . L also vanishes if ZkZi = 1 for some / 7^ k. Hence all 
the poles of L/R are removable, and since L/R is elliptic, it must be constant. To show the constant is 
nonzero, we notice that the functions inside the determinant are linearly independent, so the columns 
of the determinant are linearly independent. This concludes the proof. □ 

Remark 6.3. A more convoluted way to arrive at such determinantal representations (but the way that 
nevertheless suggested the formula above) would be to take the right hand side of the above formula 
and observe it appears in Corollary 5.4 of [War02j . What appear in the determinant on the left are the 
abelian interpolation functions i?;* discussed in the previous section: 



'p{az'^^;q)n^i (;) f") -tt . yj 0p{b/a,abq^'' ^''■,q)k- 



l<k,l<n ep{bzt';q)n-l V ^pi^^t-^Dr 
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The above formula in fact alfows us to compute the constant exphcitly by expanding the biorthog- 
onal functions in terms of abehan interpolation functions (only leading coefficient is of interest for the 
determinant, and it is explicitly given in |Rai06| ) . 

To simplify notation hereinafter we let 

^litoq") := Ri{toq';to : ti, ta, ts; WcP^i) 

The t superscript for these functions stands for the fact their arguments, as it will become apparent 
in the next proposition, are essentially locations of the particles at time t. Likewise the parameters 



depend on t [ti and Uj are implicit for t*, m* respectively; see (50) and (lOl). We'll also denote 



^lihq") = '^i{toq'')^s{tl\q, toil, tot2, tots, toUo,ptoUi)/ci so that 
Y,'^k{toq')-^iitoq')^5k,i 



s>0 



Thus Lemma 6.2 along with (181 and (50) yields: 



= const- det $ Li (io?""") • det *_i(io9'"°)- 

Kk.Kn Kk.Kn 



Proposition 6.5. We have 



s>0 



with w'q and w[ as in Theorem 
Proof. We observe that 



3.10 and Z ■ 



a ( t — lj.t-1 t-\.t-\ ^.t-\.t — \\ ■ 



s>0 



(51) 



Proposition 6.4. 

Prob{X{t) = {xi,...,xn)) = const ■ det $Li(io'?"'=) • , det *Li(iog"^) \{^., 

Kk.Kn Kk.Kn 



(52) 



which expresses the relation BA = 1 where A{k, I) = ^^(iog'), B{k, I) 



by definition (see (51)). We now apply the difference operator ©(ug ^,tQ 



'i'litoq'') and we know AB = 1 
^,t*^^) (corresponding to the 



Markov transition f i— > t — 1) to both sides and observe the parameters at time t are the required q shifts 
of the parameters at time t~ 1 (see ([48|). Finally on the right hand side we have a delta function which 
is acted upon by the difference operator to produce the desired result (see Proposition 4.6). □ 



Remark 6.6. In jBGRlO] and |BG09j formulas as in the above proposition involved discrete orthogonal 
polynomials (g-Racah and Hahn respectively) and were proven via the three term recurrence relation 
satisfied by these polynomials (which is an identity between hypergeometric or g-hypergeometric series) . 
Such a relation exists for biorthogonal functions as well (we refer the reader to |SZ00| for an explicit form, 
though with different notation) and can be used to prove the above proposition, but the computations 
are more involved. 
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Remark 6.7. A similar result holds if we apply the transition i H> i + 1 which corresponds to the 
operator ©(mi, ^2, is)- For that though, we have to renormalize the biorthogonal functions at either t2 
or ^3 (see (48) and (49)), so the bidiagonal matrix that will appear on the RHS will be of the above 
form conjugated by two diagonal matrices (coming from the renormalization coefficients). This is an 
artifact of our choice of coordinates (we are counting particles going up from the bottom left edge of the 
hexagon) . 



Finally, in applying Theorem 6.1 to the f — > t — 1 Markov chain X{t) we need to check that the 



transition probabilities have the required determinantal form. This is a consequence of Theorem 3.10 



Lemma 6.2 and the following computation (the proof of which is immediate from Theorem 3.10 and 
Proposition 6.5 we use the notation from 3.10 for w'q, w'i,X, Y): 



det {vt^t-iixk,yi)) ^ const ■ 

l<k,l<N 



(53) 



We thus obtain: 
Proposition 6.8. 

Prob{X{t - 1) = Y\X{t) =X) = const 



dcti<k,i<Nivt^t-iixk,yi))deti<kj<NiK-iiiQ 1^')) 

dcti<fe.,<Ar($*fc_l(t*g-0) 



(54) 



Theorem 6.9. The Markov processes 1 1— > i ± 1 discussed in Section^ meet the assumptions of Theorem 
6.1\ and are therefore determinantal. 



Proof. This follows from all the result s ga thered in this Secti on fo r the t— Markov chain with / = $ 



and g = 5' in the notation of Theorem 6.1 For i+ see Remark 



6.7 



□ 



Remark 6.10. For obtaining quantitative results about the artic boundary, one can try to look at the 
asymptotics of the diagonal of the correlation kernel of the process (which is of course the probability 
that a particle is present at that site): 



S+JV-l 



K{x,x) = ^ R\{toq''\tQ : ti,t2,t'i]UQ,pui)Rl{toq^\to : ti,t2,t3;pui,uo)x 



^x{tl\q, toti,tot2,tQt3, toUo,ptoUi)Ai{l/{puoUi)\q, toil, ^0^2, io^a, 1/ (^oMo), 1/ {phui)). 



7 Computer simulations 

In this section we present computer simulations of the exact sampling algorithm from Section [5] We 
are (with one exception) looking at 200 x 200 x 200 hexagons, and parameters are chosen so the elliptic 
measure sampled is positive throughout the range of the algorithm (recall that the algorithm starts 
with a 200 x 400 x box and increases c while decreasing 6 by 1, until it reaches the desired size - 
after 200 iterations in our case). Under each figure we list the values of the four parameters p, g, wi, W2. 
Computations and simulations are done using double precision, the S* i— >■ 5 + 1 algorithm polynomial 
algorithm described above, and a custom program written in Java and that can handle large hexagons 
(in excess of iV = 1000 particles) fast enough on modern CPUs. 
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Figure 11: p = q = 0.999999995, vi = 0.0000214, V2 = 1.00675. 400 x 400 x 400. Because q is very 

close to 1, the limit shape looks uniform (recall that q = 1 gives rise to the uniform measure). 



In Figure [TT] we observe that the sample looks like one from the uniform measure with the arctic 
ellipse theoretically predicted in |CLP98| clearly visible. 

Figures [12] and [13] exhibit a new behavior for the arctic circle: the curve seems to acquire 3 nodes 
at the 3 vertices of the hexagon seen in the pictures. To obtain these shapes, the parameters have been 
tweaked so that the elliptic weight ratio vanishes (or = oo) at the respective corners (in other words, 
the weight ratio Q is "barely positive" as described in Section 2.3). To be more precise, we have: 



q = eT-i 

2T-1 

vi^q 
V2 = 1/g- 

This fixes 3 of the 4 parameters of the measure and we have the extra degree of freedom p and 
so we obtain a 1-parameter family of trinodal arctic boundaries. All simulations are taken from the 



trigonometric positivity case {q,vi,V2 are of unit modulus - see Section 2.3). While the first arctic 



boundary looks like an equilateral "flat" triangle, the second looks like an equilateral "thin/hyperbolic" 



triangle. The change from Figure 12 to 13 is an increase in p (and indeed if we increase p further the 
triangle will get thinner and thinner, until it will degenerate into a union of the 3 coordinate axes as 
p — 1). The limit p — !■ yields the same "thinning behavior" in the real positivity case. 

Finally in Figure[T4]we exhibit a trinodal case in the top level trigonometric case p = when q, vi, V2 
are of unit modulus (in the case q and Vi are real, arctic boundary is the union of the coordinate axes 
as stated above). 



8 Appendix 

In this Appendix we show how one can assign S'3-invariant weights to the three types of rhombi (lozenges) 
that make up a tiling of a hexagon in the triangular lattice. We start with the 2x2x2 triangle (inside the 
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Figure 12: An instance of a trinodal arctic boundary, p — 0.00186743, arg g — 0.000835422, argwi = 
0.667502, argwa = -0.000835422. 




Figure 13: Another instance of a trinodal arctic boundary, p = 0.2,argq = 0.000835422, arg wi = 
0.667502, argW2 — —0.000835422. Note p is larger in this case than in the previous. 



triangular lattice) depicted in Figure 15 that contains an overlapping of the 3 types of rhonibi considered 
for our tilings. 



To each such type of rhombus we assign a label from the set {ui, U2, us} (see Figure 16 ) such that if 
the rhombi are as described overlapping inside a 2 x 2 x 2 triangle we have 

UiU2Ui = 1- 



Each Ui will eventually be a power of q times Ui (see Section 2.2). First, we can obviously shift 
any of such rhombi along the directions given by their edges, either upwards or downwards. If we shift 
the horizontal lozenge labeled upwards-right or upwards-left, the label of the new lozenges will be 
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Figure 14: Top level trigonometric p — case. As above, argg — 0.000835422, argwi = 0.667502, arg?;2 = 
-0.000835422. 




Figure 15: A 2 x 2 x 2 triangle composed of 3 overlapping lozenges of each type. 




Figure 16: The 3 types of rhombi (lozenges) and their labels. 

multiplied by q^^. If we shift it downwards-right /left, the label will get multiplied by q. Naturally, if 
we shift directly upwards, the label will be multiplied by q~^ as a composite of an upwards-right and 
and upwards- left shift. A similar rule is used for lozenges with labels U2 and u^. The process is depicted 
in Figure |17[ with the caveat that for labels ui and U2 we only show the directions in which the label 
gets multiplied by q (it gets multiplied by q^^ in the opposite two directions than the ones depicted). 
Clearly translating any lozenge along its long diagonal does not change its label. 
To a lozenge with label Ui (z = 1, 2,3) we assign the following weight: 

■u;i:(lozenge with label Ui) — 9p{ui), i — 1,2, 3. 
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Figure 17: Shifting lozenges in the triangular lattice, we shift the labels hy q or q ^ as depicted. 



where 



Ui = q"^ 



wi,M2,U3 are three complex numbers that multiply to 1 and {x,y,z) is the 3-dimensional coordinate of 
the center (intersection of the diagonals) of a lozenge. At this point we need to fix a choice of square 
roots: aA'^' a/"2 ' aA's such that y^y^-v/ua — 1. Note the 3-dimensional coordinates are only 
defined up to the diagonal action of Z. Figure 18 depicts the 3 lozenges with labels Ui (x = y = z = 0) 
in the chosen coordinate system. 

This way of assigning weights is manifestly 5'3-invariant. To recover the same probability distribution 



as in Section 2.2 (i.e., a gauge-equivalent weight for tilings) we again require that the weight of a tiling 
of a hexagon is the product of weights of lozenges inside it. To check this, one can simply check the 
weight ratio of a full 1x1x1 box to an empty 1x1x1 box (this is a gauge-invariant quantity) under 
the present assumptions and observe the result is the same as in ([T]). 



y 




X 



Figure 18: Choice of a coordinate system and the 3 special parameters (lozenges centered at the origin) 
needed. 

The 5*3 invariance can be viewed at the level of the partition function (the sum of weights of all 
tilings in a hexagon written in this gauge) as follows. We start with an a x /3 x 7 hexagon. The origin 
is at the hidden corner of the 3D box. In the canonical coordinates 



{ui 



qy+'-^-Ui,U2 



q-+'-^yu2,U3 



x-\-y-'2z 



U3) 



the 6 bounding edges have the following equations (see Figure 19 for correspondence between edges and 
L,'s): 
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Ug/ui := Li 
U2/U3, L2 
Ui/u2 := £3 

:= ^5 



Ui/u2 ■■= Lq := q^^ui/u2 

q^''u2/u3 
q-^"ui/u2 
q^^us/u, 

q^^^U2/u3 



(55) 




Lr 



z 

X, 



L3 




Figure 19: An a x /3 x 7 hexagon with canonical coordinates of the edges on the outside and edge lengths 
on the inside. 

We then have the following proposition. Throughout, the 53-invariant weight is assumed. 
Proposition 8.1. The partition function for an a x /3 x j hexagon is equal to: 

p ^ rp.,,,(g^+"+^+>, g^+°P, g^+^P, 

rp,g,,(qp, q^+P+lp) 

rp,g,g(g^~"+'^til, g^""+''Ml), q^-l^+^U2, q^-^+''U2,q^-"'+°'U3, q^-^+Pu^) ~ 
p ^ ^.^ Vp,,,MLoL2U?"'p, qiLoLiL5)^/^p, q{LoL,L2y/^p, (7(^2X3^4)^/^) . 



(qp, qiLo/L^y/^p, q{UlL,Y/^p, q{L2/L,y/^p) 



Tp,q,MLoL2L3f'\ qjl^L^L^Y'^ q{L2LiLr^)^/\ q{L^L2L^Y'\ qiUL^L^Y'^qjL^L^LiY'^) 



r,,,,MLo/Uy/^ q{L3/L^Y/\ qiL^/L^y/^ q{L2/Loy/^ q{U/L2Y/^ q{Li/L^Yl^) 



where 



P 



a/37- 



g ' * ^^1 "2 "3 

rp,9.t= n (i-p^+'g^+'i'+V^)(i-pVt'^)- 

j j\fc>0 
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It is left invariant by S3 permuting the coordinates Ui; equivalently, the tuple ((x, a, ui), (y, /3, U2), {z, 7, U3)). 
Furthermore, this invariance can be expanded to the group W{G2) = S3 x ^2 = Dih^ (the symmetry 
group of a regular hexagon) with the missing involution being the transformation: 

q^ui q^U2 q^U3 
where A = -2a + /3 + 7, B^a - 2/3 + 7, C = a + /3 - 27. 

Proof. We start with the elhptic MacMahon identity derived in the Appendix of |BGR10) : 

EtiiingsT^^(r^ G) ^ ^„/37 -Q epiq^+y+'-\qy+'-^-^ui,q-+'-y-^u2,q'=+y-'-^u3) 



wt{o,G) ^ 0p{q^+y+^-^,qy+^-^ui,q^+--yu2,q''+y-^u3) 

where denotes the empty tiling (box) and G is any gauge equivalent to the ones used in this paper 
(that is to say, both sides are gauge-independent). For G the 53 invariant gauge herein discussed, the 
formula for the empty tiling multiplied by the right hand side above simplifies the partition function via 
straightforward computations. We arrive at the desired result using the following transformations for T 
functions: 

rp,g,t(te) = Tp^q{x)Tp^q^tix) 

The limit p 1 is needed for technical reasons to avoid zeros of triple T functions. 

For the 5'3-invariance, it suffices to show how the edges transform under the 3-cycle (ui,U2,U3) — ^ 
(u2,M3,C(i) (a 120° clockwise rotation) and the transposition ui o U2 (a reflection in the z axis). For 
the 3-cycle, the new edges (denoted with primes) have equations: 

L'i = Lt+2 

where -t-2 is taken modulo 6, while for the transposition we have: 

L'^ = 1/^3, M = 1/^2,4 - 1/^1,4 - = 1/^5,4 - 1^4. 

Both these transformations leave the partition function invariant. The extra involution giving the 
group W{G2) is a reflection through the centroid of the hexagon having coordinates: 

{q^l^u,,q^l^U2,q^'^U3) 

The edges transform as: 

L[ = l/U+3 

where addition is mod 6. We look at the first form of the partition function written in the statement. 
We use the following two difference equations to simplify the calculations and arrive at the original form: 

^P,q,q{^l^) = ^p.q.qiPQx) = T g^g{qx)rp,g^q{qx) 



Tg^g{q q^X,x) ^ , yni -(/(",-)+m(^)) 

Tg^giq^x,q^-x) ^ 



□ 
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